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FOREWORD 

The  research  described  in  this  report  was  performed  during  the  2000-2001  fiscal  years 
as  part  of  an  effort  supported  by  the  Office  of  Naval  Research  funding  document  numbers 
N0001401WR20168  and  N0001400WR20063  to  improve  radar  classification  and 
identification  capabilities  for  noncooperative  airborne  targets.  This  continues  to  be  a 
primary  goal  of  radar  research  programs  and  considerable  effort  has  been  expended  within 
the  last  few  decades.  The  accuracy  of  any  radar  imaging  method  depends  in  a  sensitive 
way  upon  available  radar  resolution  and  signal-to-noise  ratio. 

The  current  work  describes  an  algorithm  for  determining  the  position  and  strength  of 
radar  target  scattering  centers  from  low  resolution  and  noise  corrupted  data.  The  technique 
is  remarkably  stable  against  very  poor  signal-to-noise  ratios  and  offers  good  super 
resolution  capabilities. 
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INTRODUCTION 


In  the  past  few  decades,  considerable  attention  has  been  paid  to  the  problem  of  obtaining 
accurate  and  high  resolution  radar  images  of  complex  targets  (cf,  References  1  through  3,  and 
references  cited  therein).  When  the  ultimate  goal  is  radar-based  target  recognition,  however,  even 
high  quality  radar  images  generally  contain  enough  erroneous  and  surplus  information  to  prevent 
the  use  of  practical  template-based  matching  algorithms  (References  4  and  5).  In  this  situation, 
the  conventional  heuristic  remedy  truncates  these  images  into  a  small  set  of  basic  image  “features,” 
including  points,  plates,  dihedrals,  single  and  doubly  curved  surfaces,  edges,  and  corners  (References 
2  and  4  through  7).  In  noisy  and  data-limited  environments,  practical  considerations  preclude 
distinguishing  these  scattering  mechanisms  (References  1,  6,  and  8),  and  the  most  common  set  of 
image  features  is  formed  from  estimates  of  the  position  and  strength  of  target  scattering  centers 
(assumed  to  behave  as  noninteracting  point  subtaxgets). 

But,  when  the  radar  image  is  of  low  resolution  and  contaminated  with  noise-induced  artifacts,  it 
can  also  be  very  difficult  to  accurately  estimate  either  the  position  or  strength  of  the  scattering 
centers.  Of  course,  this  estimation  error  is  directly  related  to  the  efficacy  of  template-based 
target  identification  methods  and  can  preclude  such  systems  from  application  in  many  real-world 
environments  for  which  radar  is  the  only  otherwise  viable  sensor. 

Below,  we  will  briefly  review  the  problem  of  radar  image  formation.  This  review  will  serve  to 
define  the  problem  and  establish  notational  conventions.  Using  this  framework,  we  will  then  describe 
a  method  for  accurately  estimating  the  position  and  strength  of  target  scattering  centers  directly 
from  measured  radar  data.  Examples  of  dependence  on  noise,  resolution,  and  scatterer  number 
assumptions  are  also  presented. 


RADAR  DATA  AND  THE  LINEAR  SCATTERING  MODEL 


For  notational  simplicity,  we  consider  a  monostatic  scattering  configuration  in  which  the  radar 
transmitter  is  co-located  with  the  receiver.  Denote  by  5^(0  signal  received  by  the  radar  at  time 
t  and  let  SgcattCO  response  signal  of  the  target  to  an  incident  interrogating  signal  The 

original  radar  signal  processing  problem  is  that  of  optimal  detection  in  (additive)  noise  so  that  the 
received  data  are  of  the  form:  s^{t)  =  H-  n(t),  where  n{t)  is  a  random  noise  process. 

In  conventional  radar  systems,  “optimal  detection”  has  been  accomplished  through  maximum 
likelihood  methods  that  compare  Sj.(t)  to  an  idealized  family  of  signals  predicted  using  a  model  Af. 
We  assume  that  M  is  unique  and  let  denote  the  a  posteriori  conditional  probability 

density  of  given  s^{t)  (the  a  posteriori  density  is  determined  by  the  statistics  of  n{t)  and  a 

priori  target  estimates  (Reference  9)),  Then  the  maximum  likelihood  model  satisfies 


^ML  =  arg 


max 

u€  model  space 


(1) 


It  is  usual  to  parameterize  the  model  space  to  facilitate  the  search  for  AIj^l*  active  radar, 
the  natural  parameterization  is  based  on  the  scattering  interaction  between  the  interrogating  field 
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and  the  target’s  reflectivity  behavior.  The  simplest  version  of  this  interaction  is  the  linear  radar 
(weak-scatterer)  scattering  model  that  relates  to  by 


■^scatt 


JR2 


(2) 


where  g^{t,  u)  is  the  target  reflectivity  density  defined  in  such  a  way  that  Q^{t,  u)dtdu  is  proportional 
to  the  field  reflected  from  the  target  at  range  between  ct/2  and  c(^+d^)/2  with  Doppler  shift  between 
1/  and  u  +  du.  In  general,  the  reflectivity  density  actually  depends  upon  the  incident  field.  But  when 
these  nonlinear  effects  can  be  neglected,  the  ML  model  also  yields  the  maximum  detection  signal- 
to-noise  ratio  and,  in  the  Gaussian  white  noise  and  no  prior  target  information  case,  correlation 
receivers  (References  9  through  11)  attempt  to  find  the  t  and  i/  that  maximize  the  real  part  of 

r,it,  u)=  f  s^t')  sUt'  -  (3) 


To  understand  just  what  it  is  that  the  correlation  receiver  “sees,”  substitute  the  model  of 
Equation  2  into  Equation  3  to  yield 


v(t,  ty)=  [  u')  x(t  u')  e'^''+‘''^^^-^''>^^dt'du' 

7r2 


(4) 


where 

Jr 


(5) 


and  we  have  suppressed  an  additive  noise  function  (for  the  present). 

Equation  4  is  an  imaging  equation  in  which  Q^{t,v)  is  the  object  function  and  the  remaining 
factors  in  the  integrand  represent  an  imaging  kernel  (a  point-spread  function).  Consequently, 
is  an  estimate  of  ^^(t,  i/)  and  we  can  write  ^^(t,  i/)  =  r){t^  v).  The  ideal  imaging  kernel  would  be  given 
by  x(t,  =  ^(t)^(i>'),  but  there  is  a  well-known  limit  to  the  “narrowness”  of  the  peak  of  \x{t^  y)\  that 
can  be  (simultaneously)  attained  in  the  t  and  u  directions.  This  ambiguity  relation  means  that  an 
interrogating  signal  s^^^{t)  that  leads  to  good  resolution  in  the  t-direction  will  generally  have  poor 
resolution  in  the  i/-direction,  and  vice  versa. 

A  method  for  getting  around  this  fundamental  limitation  selects  a  family  of  signals  ^^(t), 
parameterized  by  0,  which  lead  to  x${^^  ^)  f^hat  are  independently  localized  along  different  directions 
with  respect  to  target  orientation.  In  Inverse  Synthetic  Aperture  Radar  (ISAR)  imaging,  for  example, 
the  parameter  6  is  target  angular  aspect,  which  varies  in  time  if  the  target  is  rotating.  The  radar  is 
(usually)  configured  to  transmit  a  series  of  pulses  s{t—jAt)^  j  —  1, . . . ,  0,  with  good  range  resolution 
but  poor  cross-range  (Doppler)  resolution.  Assume,  for  simplicity,  a  constant  target  rotation  rate 
0  =  fi.  Assume  also  that  the  scattering  components  that  make  up  the  target  are  persistent  (i.e.. 
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independent  of  6  over  the  aperture  AO  =  {9  -  l)AtO).  Then  the  jth  pulse  will  interrogate  the 
rotated  object  function  ggit,  u)  =  p(tcos0j  +  a~^i'sin6j,  —at  sin 0^-  +i/cos9j)  where  6^  =  jAf  fl  and 
a  a  fl  is  a  scale  factor  relating  target  dimensions  in  t  to  those  in  v.  Equation  4  becomes 

r]g.{t,v)=  /  p{t'  cos 0  ■  +a“ ^ u'  sin  Oj ,  -at'  sin  0^  +  v'  cos  6^ ) 

^  Vr2  (6) 

X  I/')  dt'du’ 


Substituting  the  idealization  x(ij  t')  =  ^{t)  yields 

V9At,i')=  f  pit',v')dl  (7) 

where  dl  is  the  differential  arc  length  along  the  line  L{t;6j)  =  |  t'cos0^  +z/'sin0j*  =  t). 

Equation  7  is  the  idealized  range  profile  of  the  target  at  aspect  9  and  it  is  easy  to  see  the  relationship 
between  ISAR  imagery  and  tomographic  reconstruction  from  this  expression  (note  that  Equation  7 
is  actually  independent  of  u). 

Realizable  measurement  systems  will  obtain  finite  and  band-limited  data,  of  course,  and  will  be 
unable  to  generate  true  ^-functions  in  the  time-domain.  This  situation  has  an  idealization  that  uses 
incident  radar  signals  of  the  (frequency-domain)  form 

S{w)  =  ^{s}(a;)  =  rect  w  €  (wj,  Wj)  (8) 


where  ^  denotes  the  Fourier  transform  operation,  Au)  =  ^2  —  Wj,  Wq  ~  rect(x)  =  1 

if  X  €  (-5,5)  and  0  otherwise.  The  x{f,  i')  associated  with  this  signal  is 


X(t,  v)  =  sine  [i(Aw  -  |j/l)t]  ,  |i/l  <  Aw 


(9) 


where  sinc(x)  =  sinx/x. 

By  increasing  the  bandwidth  Aw,  the  kernel  x(<.  i')  given  by  Equation  9  can  be  made  to  more 
closely  approximate  S{t).  The  cross-range  resolution,  which  is  determined  by  AO,  cannot  generally 
be  specified  as  part  of  the  radar  system  design  and  ISAR  imaging  systems  must  often  spend 
impractically  long  time  intervals  waiting  for  the  target  to  rotate.  More  typically,  ISAR  imaging 
practitioners  simply  resign  themselves  to  low  cross-range  image  resolutions. 
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TARGET  MODEL  BASED  SIMPLIFICATION 


Considerable  further  simplification  occurs  when  the  object  function  can  also  be  replaced  by  a 
lower  dimension  parameterization.  Heuristic  arguments  (References  1  through  8)  have  been  made 
for  replacing  the  weak  scatterer  reflectivity  density  by  the  scattering  behavior  associated  with  a 
collection  of  non-interacting  characteristic  shapes.  The  simplest,  and  most  common,  choice  for 
these  shapes  is  point  scatterers  (or,  more  correctly,  point  scatterer  behavior  over  the  range  of 
frequencies  and  aspects  of  the  observation).  In  limited  and  noisy  data  situations,  the  point  scatterer 
model  is  a  practical  expedient — and  may  be  a  necessity  (Reference  6  and  8).  Assuming  this  model 
approximation  is  accurate,  we  can  write 


N 

=  “  yn) 

n-l 


(10) 


where  G  Cis  the  local  scatterer  strength  and  the  sum  is  over  all  N  scattering  centers. 
Substituting  the  model  10  into  Equation  6  (with  the  rotated  object  function)  yields 


N 


n— 1 


(11) 


where  u^{9)  =  -i„sin0  +  j/„cos0,  and  v^iO)  =  x„cos0  +  y„smd.  To  simplify  our  notation,  we 
consider  s^,  and  n  to  be  members  of  L^(T)  where  T  C  K.  Let  (x]„  =  x„,  [y]„  =  t/^,  and 

define  the  matrix 


X{t-u^,u-v„)e 


'lil'+Vn){t-Ur^)/2 


(12) 


If  we  construct  the  concatenated  matrices 


x(x,y)  = 


■xe,(x,y)' 

Xej(x,y) 

,  and  Tj  = 

.x:ee(x-y)- 

-Ve^. 

(13) 


where  [»7g](  =  i^),  then  Equation  11  can  be  written  as 


^  =  X(x,y)a  +  n 


(14) 


where  [a]„  =  a„  and  we  have  re-introduced  the  noise  function  term  for  completeness. 
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An  effective  method  for  estimating  the  model  parameters  {a„,a;„,2/n}n=i,..„iv  is  to  separate  the 
linear  components  a„  from  the  (Reference  12).  Formally,  we  can  write  the  least  squares 

solution  to  Equation  14  as 


a=  (xi(x,y)x(x,y))  ‘x^(x,y)»7  =  X  (x,y)i7 


(15) 


where  =  (x*)^  and  x“(x,y)  =  (x^(x,y)x(x,y))  ^xHx,y)  is  the  pseudoinverse  of  x(x,y). 
Consequently,  the  a  can  be  eliminated  and  an  estimate  of  x  and  y  can  be  determined  by  minimizing 


N-X(x,y)a||2  =  ||T7-x(x,y)x  (x,y)T7||2  =  ||(l-P(x,y))T7||2  (16) 


where 


P(x,y)  =  x(x,y)x  (x,y) 


(17) 


is  the  orthogonal  projection  on  the  linear  subspace  spanned  by  the  columns  of  the  matrix  x(x,y). 

Once  the  scattering  center  position  estimates  x,  y  are  obtained,  the  scatterer  strength  can  then 
be  estimated  using  Equation  15, 


NUMERICAL  ALGORITHM 


When  the  image  is  to  be  used  for  identification  purposes,  the  dimension  N  of  model  parameters 
will  usually  be  selected  to  be  small  (in  comparison  with  the  dimension  M  of  the  data  collected). 
The  problem  of  finding  the  minimum  of  Equation  16  can  be  efficiently  handled  in  terms  of  the  QR 
decomposition  (Reference  13). 

By  construction,  the  rank  of  x  ^  ^MxN  |g  case,  we  can  write 


x  =  Q 


R 

0 


(18) 


where  Q  €  is  unitary,  R  G  is  a  nonsingular  upper  triangular  matrix,  and  the  lower 

partition  is  the  (M  —  N)  x  N  zero  matrix. 

In  terms  of  this  decomposition,  the  pseudoinverse  of  x  can  be  expressed  as 


X-  =  [R-M0]Qt 


(19) 
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Substituting  the  decomposition  of  Equation  18  into  17,  it  is  easy  to  show  that 


P  =  Q 


In  0 
0  0 


Q*  and  I  -  P  =  Q 


0  0 


0  I 


M’-N  ] 


(20) 


where  is  the  N  x  N  identity  matrix.  Because  Q  is  unitary,  then  ||Qw||  =  |lw||  for  all  vectors  w 
and,  consequently,  from  Equations  16  and  20  we  have 


x,y  =  argmin 

x,y 


0  0 


0  I 


M-N 


(21) 


The  gradient  of  r(x,y)  =  ||r(x,y)||^  =  ||(I  -  P(x,y))77|p  is  given  by  (Reference  12) 


Vr(x,y)  =  -2Re{T?t(I-P)Dx(x,y)x  (^,y)v} 


(22) 


where 


(23) 


and  Xij  =  [x]ij-  terms  of  Equations  18  through  20,  we  can  write 


Vr(x,y)  =  — 2Re 


0  0 


QS)  Dx(x,y)  [R 


(24) 


If  we  denote  the  number  of  time  samples  per  aspect  by  T  (so  that  M  =  T  x  0)  then  from 
Equations  12  and  13  we  have 


-  5^*^'  +  2^j)  xit',  i'')^ 

X  sin0  +  y  cose) 

X  e‘<‘''+2*'i)‘72  cose  +  y  sine) 


(25) 


t  — +  mod  T  — lnt(t/'^)) 

^=^1+  »nt(»/T) 
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where  int(w)  denotes  the  integer  part  of  w  and  x  and  y  are  the  unit  direction  vectors. 

The  gradient  of  r  (Equation  22)  can  be  used  in  standard  iterative  schemes  for  finding  arg  min,^  y  r . 
Quasi-Newton  methods,  for  example,  set  the  j  + 1  estimate  in  terms  of  the  previous  estimate 

z^^^  by  (cf,  Reference  14) 

2(i+i)  2.^3)  ^  (26) 


where  =  H“^(z^^^)  is  the  inverse  of  the  Hessian,  =  Vr(z^-^^),  and  >  0  defines  a 

step  size.  Usually,  the  Hessian  is  estimated  numerically  and  is  chosen  so  that  r  achieves  a  local 
minimum  along  the  direction  from  z^^\ 

An  eflficient  implementation  of  Equations  16,  24  and  25,  using  Householder  transformations, 
is  presented  in  the  Appendix.  Most  of  the  elements  of  the  tensor  Dx  will  have  zero  value  and 
the  routine  accounts  for  this  by  tracking  only  the  nonzero  terms:  it  has  memory  requirements  of 
approximately  (SAT-f  1)  x  M  complex  values  plus  an  additional  2M  complex  workspace  requirement. 
The  algorithm-specific  computational  demands  are  dominated  by  the  QR  decomposition,  which 
requires  0{N‘^{M  —  N/3))  fioating  point  operations.  Of  course,  the  calculation  of  Dx  generally 
depends  on  the  form  of  x  ^md  the  implementation  (and  approximations)  that  best  match  its 
evaluation.  (In  our  case,  the  straightforward  calculation  of  Vx  that  we  employed  (cf  the  Appendix) 
accounts  for  more  than  98%  of  the  CPU  usage — implementations  more  appropriate  to  real  time 
environments  are  straightforward  (e.g.,  lookup  tables).) 


SAMPLE  RESULTS 


We  generated  synthetic  data  from  Equation  14  using  a  point  target  model  characterized 
by  X,  y,  and  a.  The  additive  noise  term  was  fabricated  using  a  Gaussian  noise  generator  and  had 
energy  determined  by  a  signal-to-noise  ratio  (SNR)  defined  so  that 


SMR  (x(x,y)a)^X(x,y)a 
ntn 


(27) 


All  of  our  examples  were  formed  firom  data  constructed  with  the  scaled  parameter  values:  i/  =  0, 
t  €  (—.5,  .5),  9  €  (— 4^,4^^),  Au  =  47r,  and  uJq  =  lOAu;  (e.g.,  roughly  consistent  with  a  3-meter  sized 
target  interrogated  with  a  signal  bandwidth  equal  to  10%  of  a  center  frequency  equal  to  2  GHz). 
The  radar  ambiguity  function  was  modeled  by  Equation  9. 

A  variety  of  nonhneax  optimization  approaches  were  tried  and  all  yielded  similar  results.  The 
figures  below  were  created  from  the  solutions  found  using  the  iteration  scheme  of  Equation  26  with 
a  Broyden-Fletcher-Goldfarb-Shanno  (BEGS)  Hessian  updating  scheme  and  an  Armijo’s  rule  line 
search  for  determining  (Reference  14). 

Figure  1  displays  the  estimates  for  four  different  values  of  SNR  (SNR  =  10,1,0.1,0.01).  The 
baseplane  {x-y  plane)  of  each  figure  contains  the  (gray-scale  encoded)  intensity  values  of  the 
traditional  (ISAR)  image,  /isar(^>2/)»  formed  by  applying  the  standard  convolution-backprojection 
algorithm  to  the  t)snR'  peak  values  of  /jsar  plotted  against  the  back  walls  of  the  figure: 
P^{y)  =  max^lJisAR(^,2/)l  is  displayed  in  the  x  =  0.5  plane;  and  Py{x)  =  max^  |/isAR(x,y)l  is 
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FIGURE  1.  A  Simple  Example  of  Scattering  Center  Estimation  in  Noisy,  Low 
Resolution  Situations  Based  on  Simulated  Radar  Data.  The  traditional  ISAR  image 
is  displayed  in  the  baseplane,  while  the  projected  ISAR  image  intensity  is  displayed 
in  the  backplanes.  The  actual  scattering  centers  are  marked  with  the  symbol  “O*” 
while  the  estimates  values  are  plotted  using  “x,”  Iteration  starting  points  are  marked 
with 


displayed  in  the  y  =  0.5  plane.  The  true  values  of  |a^|}  are  plotted  using  the  symbol  “0»” 

while  the  estimated  values  are  plotted  using  “x.”  The  seed  values  for  x  and  y  are  displayed  using 
in  the  base  plane. 

The  results  shown  in  Figure  1  are  typical  of  all  the  tests  we  performed  and  reflect  a  graceful 
localization  failure  with  increasing  noise  contamination.  Our  experience  is  that  noise  contamination 
generally  does  not  significantly  alter  the  character  of  the  hypersurface  r  but,  instead,  only  moves  the 
local  minimum  (of  course,  it  also  introduces  a  constant  offset).  Figure  2  illustrates  this  behavior  and 
is  a  plot  of  a  one-dimensional  cut  through  the  surface  r  for  different  values  of  SNR.  (The  strength 
and  location  of  the  scattering  centers  in  Figure  2  are  the  same  as  those  used  in  Figure  1.)  Because  of 
this  graceful  localization  failure  with  increasing  noise,  the  practical  limitation  of  the  method  results 
from  the  effects  of  noise  on  the  traditional  ISAR  image,  because  a  low  quality  image  can  prevent  an 
estimate  of  N  and  the  seed  values  for  x  and  y.  (In  general,  we  found  the  algorithm  to  be  robust 
against  seed  value  choice,  although  solution  speed  was  somewhat  altered.) 
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FIGURE  2.  The  Cost  Function  for  Various  SNRs.  Plotted  are  1-D  cuts  through  r. 


When  the  SNR  is  favorable,  however,  the  method  can  be  an  effective  localizer  and  will  usually 
sinpass  the  standard  resolution  limitations  of  an  ISAR  image.  Figure  3(a)  shows  the  “super¬ 
resolution”  capabilities  of  the  algorithm.  The  x  and  y  resolution  cells  associated  with  these  data 
are  Ax  w  2wc{u>qA6)~^  w  0,36  and  Ay  w  nc(2Au;)~^  «  0.12,  respectively.  The  actual  scatterers 
were  separated  by  Xj  -  X2  =  0.05  «  Ax/7  and  yi  -  y2  =  0.03  w  Ay/4.  Of  course,  super-resolution 
requires  that  the  user  have  some  prior  reeison  to  suspect  that  there  is  more  than  one  scatterer  in  a 
particular  target  region — this  information  surely  cannot  be  divined  from  the  raw  ISAR  image,  and 
the  model  order  N  must  be  set  in  advance  (along  with  seed  guesses).  Figure  3(b)  shows  the  results 
of  estimating  N  to  be  larger  than  its  correct  value  (in  this  case,  =  2,  =  1).  Observe  that 

the  algorithm  correctly  fits  one  of  the  scatterers  to  the  data  and  assigns  a  strength  of  lojl  =  0  to 
the  other.  (It  is  possible  for  the  method  to  converge  to  different  results,  but  the  behavior  of  Figure 
3(b)  was  the  most  commonly  observed.) 


CONCLUSION 


Traditionally,  even  the  simple  problem  of  fitting  a  point  scatterer  model  to  ISAR  data  has 
frequently  proven  impracticable  when  time  and  computer  assets  are  Umited.  By  taking  advantage  of 
the  separabUity  of  the  problem,  we  have  shown  how  the  computational  complexity  can  be  significantly 
reduced.  Moreover,  we  have  constructed  a  memory-efficient  implementation  and  demonstrated  it 
using  synthetic  data. 

The  method  is  robust  agmnst  noise  contamination  and  displays  good  super-resolution 
capabilities.  The  number  of  scattering  centers  to  be  localized  is  a  user  defined  input  (along  with 
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(6) 


FIGURE  3.  Super-Resolved  Scattering  Centers.  In  Figure  3(a)  the  point  scatterers 
are  separated  by  approximately  0.14Ax  and  0.25Ay  (where  the  resolution  cell  is 
Ax  X  Ay  in  size).  Figure  3(b)  shows  the  results  of  inappropriately  associating  two 
scatterers  to  a  location  where  only  one  true  scatterer  actually  exists. 


seed  estimates  based  on  low  resolution  ISAR  images)  and  the  “best  fit”  location  and  strength  of  the 
scatterers  are  estimated  according  to  Equations  15  and  21.  Consequently,  the  approach  appears  to 
be  well-suited  to  radar-based  target  identification  problems  that  rely  on  template  matching  schemes 
in  which  only  the  information  associated  with  the  N  strongest  scatterers  are  used. 

Throughout,  we  have  concentrated  on  the  point  scatterer  model,  but  the  method  can  be  easily 
modified  to  include  more  complex  cases  (References  6  through  8).  In  many  practical  situations, 
however,  data  limitations  and  noise  contamination  issues  can  be  expected  to  make  such  modifications 
ineffectual. 
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Appendix 

COMPUTER  CODE  FOR  CALCULATING  r  AND  Vr 


The  following  computer  code  was  implemented  in  MATLAB  and  evaluates  Equations  16,  24, 
and  25  for  x  defined  by  Equation  9.  The  algorithm  is: 


1.  Find  X  =  QJi  and  set  r  =  +  1  :  M)|p. 


2. 


Calculate  x  =  [R  ^  1 0]^  and  C  —  Q 


0  0 

.0  Im-n 


i). 


3. 

4.  Vr  =  — 2Re(4fa;). 


function  [r,  Vr]  =  rjir{z) 
global  7)  N  M  T  0  Acj  ljq  a 

^  =  [01 ,  ^2*  •  ■  -j^  vector  of  angles 
T  =  number  of  time  measurements  at  each  Oj 
N  =  number  of  scattering  centers 
M  =  number  of  observations  (=  T  x  length (0)) 
z  =  [xi,X2,  ...,a:;\^,2/l>2/2,  — 

Vx  =  M  X  2N  array  of  derivatives 

Note,  only  the  nonzero  components  are  stored. 

X  =  zeros(M,  N);  Vx  =  zeros(M,  2N); 

V  =  zeros(M,  1);  =  zeros(Ar,  1);  x  =  zeros(N,  1);  Vr  =  zeros(2iV,  1); 

I  Obtain  the  orthogonal  factorization  of  x  ^  X  =  ^  and  4^  =  <?♦  TJ. 

[x.  Vx]  =  chiJ3chi(z,  x,  Vx);  V’  =  »?; 
for  I  =  1  :  N 

[v(i  :  M),/3(i)]  =  house{x{i  •  A/,  i)); 

x(i  :  M^i  :  N)  =  row-house{x{i  •  ‘  N),v{i  : 

:  M)  =  rowJiouse{%l>{i  :  A/),u(i  :  M),y9(i)); 

I  Save  the  Householder  vectors  in  the  lower  triangle  of  x  (for  later  use) 

X{i  +  1  :  M,  i)  =  v(t  -I- 1  :  M); 
end 

I  Calculate  the  residual, 

r=|WN  +  l:M)||2; 

Computation  of  a  (by  Backsubstitution): 
a  =  [R{1  :  TV,  1  :  AT)-!  ^’(l  :  Ar);zeros(M  -  TV,  1)]; 

(a(l  :  TV)  is  also  the  scatterer  strengths.) 

a(TV)=V'(TV)/x(TV,TV); 
for  i  =  TV  —  1  :  — 1  :  1 

a(i)  =  (7/>(i)  -  x(i,  i  +  1  :  TV)  a(z  +  1  :  N))/x{h  i)\ 
end 

C  =  (/  —  P)t]  =  Q  [zeros(TV,  l);i/»(TV  +  1  :  M)]. 

Note,  Q  is  stored  in  factored  form:  Q  =  Qi  x  (32  x  Qs  •  •  •  Qn» 
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I  SO  ~  . . .  Q3  X  Q2  X  Q\. 

C  —  [zeros(Ar,  1);  +  1  :  A/)]; 

for  i  =  AT  :  -1  :  1 

v(i)  =  1;  v(i  4- 1  :  M)  =  xi^  -f  1  :  M,  i); 

C(i  :  M)  =  rowJiouse(C(i  :  M),v(i :  M),conj(/?(i))); 
end 

I  ^  =  C^Vx 

=  zeros(2N,  AT); 
for  i  =  1  :  AT 

i)  =  Ct  Vx(:,  i);  '*'(»  +  N,  t)  =  Ct  Vx(:,  i  +  Af); 

end 

I  Finally,  put  it  all  together  to  calculate  the  gradient. 

Vr  =  “2real(^a); 

function  [v,  0\  =  house{x) 

Compute  the  Householder  vector. 

(cf.  §5.2.10  of  Golub  and  Van  Loan,  third  ed.) 

s  =  ||x|{  (sign(i(l))  +  (x(l)  ==  0));  |  require  sign(O)  =  1. 

V  =  x\ 

if  s  ==  0,  /?  =  1;  return,  end 
v(l)  =  t;(l)  +s; 

function  A  =  rawJiouse{Ay  v,  /?) 

I  Overwrites  A  with  QvA  where  Qv  =  I  —  • 

w  =  A^  v\ 

A  =  A^  pvw^\ 

function  [x,  Vx]  =  chiJDchi{z,  x,  Vx) 
global  Tj  N  M  T  0  Alj  loq 

I  Calculate  x  Vx- 

ct  =  cos(0);  st  =  sin(0); 

(/  =  -2(1  :  AT)  st^  +  z{N  +  l:  2N)  ct^;  V  =  z{l  :  N)ct^ z{N  +  1  :  2N)  st^] 

for  fc  =  1  :  AT 
for  j  =  1  :  M 
je  =  ceil(j/r); 

t  =  (1  +  mod{j  -  l,r))/T  -  .5  -  U(kJe); 

[x»  ^x]  =  oTn6-Z>om6(t,  u); 

e  =  exp(— ii/t/2);  |  note:  i/  +  2v  =  v  (=  — i/) 

=  xe; 

a=  -L>x(l)  +  ii"X/2; 

6=-Z>x(2)  +  itx/2; 

Vx(i,fc)  =  (-ast(jfl)  4-6ct(j©))e; 

Vx0‘,  N  *f  fe)  =  -(act(jfl)  4-  6st(jfl))e; 
end 
end 

function  [x,  Dx]  =  om6_Dam6(t,  i/) 
global  T)  N  M  T  0  Aw  wo 

I  Calculate  x»  ^x/^f  and  dx/du 

6  —  Aw  —  1j/|;  spni,  =  sign(i/)  4-  (i'  ==  0); 
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if  <5  <  0,  X  =  0;  =  0;  return,  end 

arg  —  6 1/2;  e  =  exp(i(a;o  + 

C  =  e^/Au;;  sine  —  smc(arp); 

I  Determine  the  ambiguity  function 

X  =  C  sine; 

I  Determine  the  derivatives 
if  \arg\  >  0 

sind  =  (cos  (arp)  —  sinc)/arg\ 
else 

sine'  =  0; 
end 

if  (5  >0 

Dx(l)  =  i(ct^o  +  i^/2)x  +  C(5sine'/2; 

Dx(2)  =  (it/2  —  sgnulS)x  -  Cs^n^  t  sme'/2; 
else 

I>x(l)-0; 

Dx(2)  =  -espn,y/Aa;; 
end 
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